#include "stdafx.h"
#include "PNGFinderSimulationAlgorithm.h"

using namespace System;
using namespace System::Collections::Generic;
using namespace System::Linq;

namespace SpikingNeuronLib
{
	PNGFinderSimulationAlgorithm::PNGFinderSimulationAlgorithm(List<Connection^>^ connections, int numExcitatoryNeurons, int totalNeurons, int groupSize, int maxPNGLength)
		: PNGFinder(connections, groupSize)
	{
		_a = gcnew array<double>(totalNeurons);
		for (int i = 0; i < numExcitatoryNeurons; i++)
			_a[i] = 0.02;

		for (int i = numExcitatoryNeurons; i < totalNeurons; i++)
			_a[i] = 0.1;

		_d = gcnew array<double>(totalNeurons);
		for (int i = 0; i < numExcitatoryNeurons; i++)
			_d[i] = 8;

		for (int i = numExcitatoryNeurons; i < totalNeurons; i++)
			_d[i] = 2;

		_virtualEvents = gcnew array<VirtualEvent, 2>(maxPNGLength, totalNeurons);
	}

	int PNGFinderSimulationAlgorithm::SearchCurrentTriplet(String^ tag, int anchorNeuron, array<int>^ triplet, int maxDelay, int numPolychronousGroups,
		int groupSize, double weightLimit, int maxPNGLength, int numExcitatoryNeurons, int totalNeurons,
		Dictionary<Tuple<int, int>^, Connection^>^ connectionPairDictionary,
		Dictionary<int, Dictionary<int, List<Connection^>^>^>^ groupedConnectionDictionary, int verboseLevel)
	{
		return SearchCurrentTripletSimulated(tag, anchorNeuron, triplet, maxDelay, numPolychronousGroups,
			groupSize, weightLimit, maxPNGLength, numExcitatoryNeurons, totalNeurons,
			connectionPairDictionary, groupedConnectionDictionary,
			NumPostspikes, GroupEvents, _virtualEvents, Links, _a, _d, verboseLevel);
	}

	int PNGFinderSimulationAlgorithm::SearchCurrentTripletSimulated(String^ tag, int anchorNeuron, array<int>^ triplet, int maxDelay, int numPolychronousGroups,
		int groupSize, double weightLimit, int maxPNGLength, int numExcitatoryNeurons, int totalNeurons,
		Dictionary<Tuple<int, int>^, Connection^>^ connectionPairDictionary,
		Dictionary<int, Dictionary<int, List<Connection^>^>^>^ groupedConnectionDictionary,
		array<int>^ numPostspikes, array<FiringEvent>^ groupEvents, array<VirtualEvent, 2>^ virtualEvents,
		List<LinkedEventPair^>^ links, array<double>^ a, array<double>^ d, int verboseLevel)
	{
		const int minPNGPathLinks = 7;

		// MG originally used unassigned variable 'timing' that was assumed to be set to the max timing for each iteration
		// based on the order of iterating through the delays
		// this new implementation cannot assume that delays are iterated in a particular order
		// therefore the new implementation explicity sets the max timing for each outer iteration of time, t
		int maxTiming = 0;
		int latency = maxDelay;			// maximum latency 
		int numFired = groupSize;

		// reinitialise
		links->Clear();

		array<int>^ layer = gcnew array<int>(maxPNGLength);
		array<double>^ v = gcnew array<double>(totalNeurons);
		array<double>^ u = gcnew array<double>(totalNeurons);
		array<double>^ I = gcnew array<double>(totalNeurons);

		// fetch the triplet connections
		array<Connection^>^ tripletAnchorConnections = gcnew array<Connection^>(triplet->Length);
		for (int i = 0; i < triplet->Length; i++)
		{
			Connection^ connection = connectionPairDictionary[gcnew Tuple<int, int>(triplet[i], anchorNeuron)];
			tripletAnchorConnections[i] = connection;
		}

		// find the longest connection for the triplet
		int Dmax = 0;
		for (int i = 0; i < groupSize; i++)
		{
			Connection^ connection = tripletAnchorConnections[i];
			if (Dmax < connection->Delay)
			{
				Dmax = connection->Delay;
			}
		}

		// initialize the data structures
		// for each triplet neuron...
		for (int groupIndex = 0; groupIndex < groupSize; groupIndex++)
		{
			int neuronIndex = triplet[groupIndex];
			Connection^ anchorConnection = tripletAnchorConnections[groupIndex];

			// add the neuron (and firing time) to the new PNG
			groupEvents[groupIndex] = FiringEvent(Dmax - anchorConnection->Delay, neuronIndex, EventLabel::Foreground);

			// the layer for each neuron is one more than the layer of its presynaptic neuron in the group (starting from layer 1 for the initial triplet)
			layer[groupIndex] = 1;

			// for each outgoing connection from neuronIndex
			Dictionary<int, List<Connection^>^>^ outgoingConnections = groupedConnectionDictionary[neuronIndex];
			for each(int delay in outgoingConnections->Keys)
			{
				for each(Connection^ connection in outgoingConnections[delay])
				{
					// if the connection is to a strong excitatory neuron that will fire after the anchor...
					if ((connection->Weight > weightLimit) && (delay >= anchorConnection->Delay))
					{
						// found a candidate that has a sufficiently weighty synapse and spikes after i
						// tentatively add to the group
						int timing = groupEvents[groupIndex].Time + delay;
						if (timing > maxTiming) maxTiming = timing;

						VirtualEvent vevent = VirtualEvent(
							groupEvents[groupIndex].NeuronIndex,    // presynaptic neuron
							connection->PostNeuron,					// postsynaptic neuron
							delay,
							connection->Weight,
							0);										// link level unused

						virtualEvents[timing, numPostspikes[timing]++] = vevent;
					}
				}
			}
		}

		// initialise membrane potential (v) and recovery vectors (u), and the input vector (I)
		for (int neuronIndex = 0; neuronIndex < totalNeurons; neuronIndex++)
		{
			v[neuronIndex] = -70;
			u[neuronIndex] = 0.2 * v[neuronIndex];
			I[neuronIndex] = 0;
		}

		// fire all neurons in the current triplet set and their downstream targets
		int t_last = maxDelay + maxDelay + latency + 1;
		int t = -1;
		while ((++t < t_last) & (numFired < maxPNGLength))		// t_last is a variable time limit
		{
			// for each candidate, add the connection weight to the postsynaptic input
			for (int p = 0; p < numPostspikes[t]; p++)
			{
				I[virtualEvents[t, p].PostNeuron] += virtualEvents[t, p].Weight;
			}

			// update v and u for all neurons and reset the input
			for (int neuronIndex = 0; neuronIndex < totalNeurons; neuronIndex++)
			{
				v[neuronIndex] += 0.5 * ((0.04 * v[neuronIndex] + 5) * v[neuronIndex] + 140 - u[neuronIndex] + I[neuronIndex]);
				v[neuronIndex] += 0.5 * ((0.04 * v[neuronIndex] + 5) * v[neuronIndex] + 140 - u[neuronIndex] + I[neuronIndex]);
				u[neuronIndex] += a[neuronIndex] * (0.2 * v[neuronIndex] - u[neuronIndex]);
				I[neuronIndex] = 0;
			}

			// find neurons that have fired and add them to the candidate list
			for (int neuronIndex = 0; neuronIndex < totalNeurons; neuronIndex++)
			{
				if (v[neuronIndex] >= 30)
				{
					v[neuronIndex] = -65;
					u[neuronIndex] += d[neuronIndex];

					if (numFired < maxPNGLength)
					{
						// update the firing event list
						groupEvents[numFired++] = FiringEvent(t, neuronIndex, EventLabel::Foreground);

						// for each outgoing connection from neuronIndex
						Dictionary<int, List<Connection^>^>^ outgoingConnections = groupedConnectionDictionary[neuronIndex];
						for each(int delay in outgoingConnections->Keys)
						{
							for each(Connection^ connection in outgoingConnections[delay])
							{
								// if the synapse is strong or the neuron is inhibitory
								if ((connection->Weight > weightLimit) || (neuronIndex >= numExcitatoryNeurons))
								{
									// found a candidate that has a sufficiently weighty synapse and spikes after i
									// tentatively add to the group
									int timing = t + delay;
									if (timing > maxTiming) maxTiming = timing;

									// BUGFIX 1: it is possible for 'timing' to exceed the maximum index of the virtualEvents arrays
									// BUGFIX 2: it is possible for 'numPostspikes[timing]' to exceed the maximum index of the virtualEvents arrays
									// the next line check these limits
									if ((timing < maxPNGLength) && (numPostspikes[timing] < totalNeurons))	// check the limit
									{
										VirtualEvent vevent = VirtualEvent(
											neuronIndex,                // presynaptic neuron
											connection->PostNeuron,		// postsynaptic neuron
											delay,
											connection->Weight,
											0);							// link level unused

										virtualEvents[timing, numPostspikes[timing]++] = vevent;
									}
								}
							}
						}

						// push back the time limit if necessary
						if (t_last < maxTiming + 1)
						{
							t_last = maxTiming + 1;
							if (t_last > maxPNGLength - maxDelay - 1)
							{
								t_last = maxPNGLength - maxDelay - 1;
							}
						}
					}
				}
			}
		}

		// if the number of fired neurons is over the threshold, calculate links and layers
		if (numFired > 2 * groupSize)
		{
			int maxLinkLevel = 0;

			// for each fired event (other than the events from the original triplet)
			for (int i = groupSize; i < numFired; i++)
			{
				layer[i] = 0;

				// find convergent virtual events that may have been causative in each firing event
				for (int p = groupEvents[i].Time; (p > groupEvents[i].Time - latency) && (p >= 0); p--)
				{
					for (int j = 0; j < numPostspikes[p]; j++)
					{
						if ((virtualEvents[p, j].PostNeuron == groupEvents[i].NeuronIndex) && (virtualEvents[p, j].PreNeuron < numExcitatoryNeurons))
						{
							// determine which layer the current event belongs to (the first layer is 1)
							for (int k = 0; k < i; k++)
							{
								if ((groupEvents[k].NeuronIndex == virtualEvents[p, j].PreNeuron) && (layer[k] + 1 > layer[i]))
									layer[i] = layer[k] + 1;
							}

							// add the new link
							VirtualEvent vevent = virtualEvents[p, j];
							LinkedEventPair^ eventPair = gcnew LinkedEventPair(
								FiringEvent(p - vevent.Delay, vevent.PreNeuron, EventLabel::Foreground),
								FiringEvent(groupEvents[i].Time, groupEvents[i].NeuronIndex, EventLabel::Foreground),
								vevent.Weight,
								vevent.Delay,
								layer[i]);

							links->Add(eventPair);

							if (maxLinkLevel < layer[i])
							{
								maxLinkLevel = layer[i];
							}
						}
					}
				}
			}

			// if the maximum layer reaches a threshold...
			if (maxLinkLevel >= minPNGPathLinks)
			{
				// discard the group if any member of the original triplet has no links to excitatory neurons
				bool discard = false;
				array<int>^ used = gcnew array<int>(groupSize);
				for (int i = 0; (discard == false) && (i < groupSize); i++)
				{
					used[i] = 0;
					for (int j = 0; j < links->Count; j++)
					{
						if ((links[j]->PreEvent.NeuronIndex == groupEvents[i].NeuronIndex) && (links[j]->PostEvent.NeuronIndex < numExcitatoryNeurons))
						{
							used[i]++;
						}
					}

					if (used[i] == 1)
					{
						discard = true;
					}
				}

				if (discard == false)
				{
					numPolychronousGroups++;

					if (verboseLevel >= 1)
					{
						Console::WriteLine("{0} {1} {2} {3} {4} {5}", tag, numPolychronousGroups, anchorNeuron, numFired, maxLinkLevel, groupEvents[numFired - 1].Time);
					}

					// throw a 'found' event or write to the console
					this->OnFoundPNG(tag, numPolychronousGroups, anchorNeuron, numFired, maxLinkLevel, groupEvents[numFired - 1].Time, groupEvents, links);
				}
			}
		}

		// reinitialise the values of N_postspikes in the utilised range to zero...
		for (int dd = Dmax; dd < t_last; dd++)
			numPostspikes[dd] = 0;

		// ...and any additional values where the time limit has been pushed back
		if (t_last == maxPNGLength - maxDelay)
		{
			for (int dd = t_last; dd < maxPNGLength; dd++)
				numPostspikes[dd] = 0;
		}

		return numPolychronousGroups;
	}
}
